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Abstract 

We study neutron matter by combining pionless effective field theory with non-perturbative 
lattice methods. The neutron contact interaction is determined by zero temperature scattering 
data. We simulate neutron matter on the lattice at temperatures 4 and 8 MeV and densities 
below one-fifth normal nuclear matter density. Our results at different lattice spacings agree with 
one another and match bubble chain calculations at low densities. The equation of state of pure 
neutron matter obtained from our simulations agrees quantitatively with variational calculations 
based on realistic potentials. 
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I. INTRODUCTION 



The equation of state of dilute neutron matter is of central importance to the structure 
and evolution of neutron stars . In addition to that, the neutron matter problem has 
many interesting physical aspects. Neutron matter has positive pressure at all densities 
and becomes a superfluid at sufficiently small temperature. If the density is small, pairing 
is expected to take place in an S-wave, but at higher density P-wave pairing might be 
dominant |3]. The neutron matter problem contains a number of very different scales. The 
neutron scattering length is very large, a„„ ~ — 18 fm, which implies that the dimensionless 
parameter kplannl ^ 1 for densities p > 10~'^pn- Here, kp = (Svr^p)^/^ is the Fermi 
momentum and ~ 0.17 fm~^ is the saturation density of nuclear matter. The effective 
range, on the other hand, is of natural size, r„„ ~ 2.8 fm. As a consequence, the parameter 
kplfrml is neither large nor small for densities p ^ p^. 

If the density is very small, p < O.lpAr, then kp\rnn\ is a small parameter and neutron 
matter is close to the limit in which kplannl —>■ oo and kplrnnl — > 0. In this limit dimensional 
analysis implies that the energy per particle and the gap have to be proportional to the Fermi 
energy 

A ^5 2m' ^2m ^ ^ 

The determination of the two dimensionless parameters ^ and C is a fascinating non- 
perturbative problem that has received a lot of attention recently. This interest was fueled 
by experimental advances in creating cold, dilute gases of fermionic atoms tuned to be near 
a Feshbach resonance jj, 0, 0, Q, • 

The traditional approach to the neutron matter problem is based on the assumption that 
nucleons can be treated as non-relativistic point-particles interacting mainly via two-body 
potentials. The two-body potentials are fitted to experimental data on nucleon-nucleon 
scattering. The many-body problem is addressed by solving the many-body Schrodinger 
equation using variational methods or Green function Monte Carlo methods guided by vari- 
ational wave functions 0, Q| . 

Even though this method has been very successful it is desirable to seek an alternative 
approach that is more directly related to QCD, systematically improvable, and that lends 
itself to numerical studies which do not rely on variational wave functions. Such an approach 
is provided by effective field theory. The use of effective field theory (EFT) methods in 



2 



nuclear physics was pioneered by Weinberg Over the last few years EFT methods 

lave been applied successfully to the study of two and three-body systems at low energy 
1^ lis! 11^ . Nuclear and neutron matter was studied using a perturbative expansion in 

powers of the Fermi momentum and using lattice simulations 



16, 



17|. 



In this paper we shall study dilute neutron matter using a nuclear effective field theory on 
the lattice. Since we are interested in densities below nuclear matter saturation density we 
shall assume that the relevant momenta are smaller than the pion mass and that we can use 
an effective field theory that contains only neutrons. In this work we shall limit ourselves to 
the lowest order effective Lagrangian which contains a single four-fermion contact interaction 
with no derivatives that is adjusted to the neutron-neutron scattering length. This effective 
theory is sufficient in order to investigate universal properties in the limit fc_F|a„„| oo, 
kplrnnl 0. An important advantage of the model is the fact that in the case of an 
attractive interaction there is no sign problem at finite density As a consequence, the 
theory can be simulated efficiently using standard hybrid Monte Carlo algorithms 

The paper is organized as follows. In Sects. II-IV we introduce the lattice theory. In 
Sect. V we discuss how to determine the coefficient of the four fermion interaction by match- 
ing to the two-body scattering length. In Sect. VI we study a low density approximation to 
the partition function based on summing particle-particle chains. In Sect. VII we describe 
our hybrid Monte Carlo method. Numerical results for the neutron density, the energy per 
particle and the equation of state are given in Sects. VIII-XII. 



II. NOTATION 



Before describing the physics we first define some notation we use throughout our discus- 
sion. We let n represent integer- valued lattice vectors on our 3-1-1 dimensional space-time 
lattice. We use a subscripted "s" such as in to represent purely spatial lattice vectors. 
We use subscripted indices such as i,j for the two spin components of the neutron, f and 
|. We let be the unit lattice vector in the time direction and let = 1, 2, 3 be the 
corresponding unit lattice vectors in the spatial directions. A summation symbol such as 

E (2) 

Is 

imphes a summation over values = 1, 2, 3. 
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We take the neutron mass to be 939 MeV, and normal nuclear matter density to be 0.17 
fm~^. We let a be the lattice spacing in the spatial direction and L be the number of lattice 
sites in each spatial direction, at is the lattice spacing in the temporal direction and Lt is 
the number of lattice sites in the temporal direction. We let at be the ratio between lattice 
spacings, 

= f ■ (3) 
Throughout we use dimensionless parameters and operators, which correspond with physical 
values multiplied by the appropriate power of a. In the end, however, we report final results 
in physical units such as MeV or fm~^. In cases where there may be confusion, we use the 
subscript phys to identify quantities in physical units. 

We use a, to represent annihilation and creation operators for the neutron, whereas 
c, c* indicate the corresponding Grassmann variables in the path integral representation. 
We let mjv be the mass of the neutron and be the neutron chemical potential. For the 
neutron fields we apply periodic boundary conditions in the spatial directions and antiperi- 
odic boundary conditions in the temporal direction. For each neutron momentum we use 
the notation 

h^Cf^ko,'fh,'fk2,'fh), (4) 

where ki, A;2, and are integers and ko is an odd half-integer. In physical units the 
momentum is 

oO^S/c^ia ^,k^2(i ^,k^3a ^). (5) 

Unless otherwise indicated, our the momentum labels will follow this convention. For 
convenience we also define 

h = (6) 

and 

= 6/i - 2/i^cos(A;*/J. (7) 

Is 

We let D^'^^^{k)5ij be the free neutron propagator. For notational convenience the spin- 
conserving Sij in the neutron propagator will be implicit. The self-energy, S(/c), is defined 
by 

Df-i\k) ^1^''^^^ ^ . (8) 

1 - j:{k)Df^-'-ik) 



where D^'^'-\k) is the fully-interacting propagator. 

In our plots we use the abbreviation "fc" for free continuum results, "f ' for free lattice 
results, "b" for bubble chain calculations, and "s" for lattice simulations results. In addi- 
tion to these abbreviations, we will use the shorthand labels shown in Table 1 for various 
combinations of spatial and temporal lattice spacings presented in our analysis. 

Table 1: Shorthand labels for various lattice spacings used 



a-i(MeV) 


a,-^(MeV) 


Label 


50 


24 





60 


32 


1 


60 


48 


2 


70 


64 


3 


80 


72 


4 



III. FREE NUCLEON 



On the lattice the free neutron Hamiltonian can be written as 
Hnn = [("^^ ~ ^ + ^)4i^s)ai{ns) 

ns,i 

~ 2^ 5Z [4i^s)ai{ns + Is) + a\{ns)ai{ns - I, 

We can approximate the partition function as a Euclidean lattice path integral, 

Z^-^ = Tr exp [-(3Hfi^] ~ zl'^' J DcDc* exp [S^"-''] , 
where Zq^"^^ is a constant and 

Sfree ^ ^ [c* {n) c,{n + 6) - e-("^^-^)"* (1 - 6h)c*{n)ci{n)] 

n,i 



nis.i 



(9) 



(10) 



(11) 



We have taken a slightly different form than that used in |17|1. Instead of the e that 

n 

appears in [17[, we use the more standard 1 — 6/i as the coefficient multiplying c*{n)ci{n). 
It is conventional to define a new normalization for q. 



-{mff-ij,)at 



(12) 



Then 



where 



Z!:^" ~ ^/'^-^e-'^"*^-'^)^^' j Dc'Dc* exp [-S^"-''] 



Sfree ^ ^ [e^'^^-'^^^'c* (n)c;(n + 0) - (1 - 6/i)c*(n)c;(n)] 



In momentum space we have 

k,i 



The free neutron correlation function on the lattice is 

/ Dc'Dc*c'i{n)c*{0)exp [S^''^^] _ 1 



/ Dc'Dc* exp [Sfr 



(no sum over i) where the free neutron propagator is 

1 



Df'^'ik) 



g-ifc,o+(m;v-M)«t - (1 - 6h) - 2hY^i^ cos{Ki^) 
1 



j-ife*o+("T'jv-/i)a;t _ ]^ _|_ 



(13) 



(14) 



(15) 



(16) 



(17) 



IV. NEUTRON CONTACT TERM 



There are two contact interactions at lowest order in the effective theory of nucleons 
without pions. But since we are considering pure neutron matter, this reduces to one 
contact interaction of the form 



Hnnnn = C^a\{ns)a'i{ns)a\{ns)ai{ns 



(18) 



Since 



exp 



Cat 



— / (is exp — -s^ + sV— Cci;(a|a-|- + ajaj^) 




(19) 



we can write 



exp 



fT r 



ds exp 



--5^+ ( sV-Ca + 



^ ) (a{«T + ^t'^i) 



(20) 
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With this interaction the partition function can be approximated by 

Zg = Tr exp + ff^^viVTv)] - j DsDcDc* exp [-S] , (21) 

where zq is a constant and 

S = Y^ [e("^^-'^)"*c*(n)c',(n + 0) - ev'^^^(")+^(l - 6/i)c*(n)c',(n) 
- 5^ [c:(n)c:(n + g + c*(n)c:(n - + 1 ^ (22) 

n,ls,i n 

This lattice action is quite simple, and in the future it may be worth considering improved 
actions in order to reduce discretization errors. Nevertheless our lattice action maintains 
some important properties. One property is that the chemical potential, /i, is coupled to 
an exactly conserved neutron number operator. This is clear since /i appears in the same 
manner as a temporal gauge link. Another feature is that in the limit as ttiat ^ oo, we find 

Tr exp [-P{Hnn + Hnnnn)] = zq J DsDcDc* exp [-S] + 0{m]^^). (23) 

Therefore any dependence on the temporal lattice spacing is suppressed by a factor of m^^. 
This makes it possible to take the static neutron limit as a precision test of the simulation 
results. We have found this test quite useful in the process of code development and 
checking. 

V. DETERMINING COEFFICIENTS 

The interaction coefficient C must be determined for various lattice spacings a and at- 
We do this by summing all bubble chain diagrams contributing to neutron-neutron scattering 
as shown in Fig. ^ 

The next step is to locate the pole in the scatte ring amplitude and compare with Liischer's 
formula for energy levels in a finite periodic box |2fl[ I21I1. 

^°-^;^^^~'^~r + '^"L^ + '"J' ^^^^ 
where ci = —2.837297, C2 = 6.375183. We then tune the coefficient C to give the physically 
measured ^5*0 scattering length. Since this scattering length is much larger than any other 



£>X) - o( 

n = 0^ 12 n \ 



FIG. 1: Bubble chain diagrams contributing to neutron- neutron scattering. 



length scale, we are in essence probing the universal behavior of interacting fermions at 

1 c 

infinite scattering length. For our results we have used a^^att — ~24 fm, though using the 
value of —18 fm specific for neutron-neutron scattering changes the operator coefficient by 
only 1%. 

The full bubble chain will have a pole when the amplitude for a single bubble times one 
vertex coefficient equals 1. We take the center of mass frame and let the total incoming 
momentum of the two neutrons in physical units be 

Pphys = (p*oa^~^ 0^ 0)- (25) 

Since the physical pole occurs in Minkowski space, in the end po will be imaginary. If we 
set = 0, then the amplitude for one bubble times one vertex coefficient is 

(l-6/i)2(e-^"*-l)i?(po), (26) 

where 

k 

and the condition for the location of the pole is 

^(^») = (l-6ft)MU°--l) - 

By the definition of Uk in ((Tj) we see that < cj^ < 12h. We assume that the lattice 
spacing in the temporal direction is sufficiently small so that h < ^. In practice this presents 
no problem since m^v is quite large. We then have 

< Wfc < 2, (29) 
-1 < 1 - cjfc < 1. (30) 
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We now make a variable transformation, 



= e . 



(31) 



We also take the zero temperature limit, — > oo, and convert from the discrete sum over 
ko to an integral clockwise over the unit circle in z using 

2n 

dz — —i—zdho, (32) 
Lt 



dkr, — i — —dz. (33) 



We then find 



, i sr^ f dz 



27rL3 ^ / ^ (e"*^"*e-*P*o/22; - 1 + a;fe) (e"*Jvate-»p.o/22;-i _ 1 + 



E 



dz 



27rL3 ^ / (e'"~"*e-'P*o/2^ - (1 - Wfe)) (e^'Jvate-ip^.o/a _ (^1 _ ^) 

27rL3 ^^^^^ J {z- e-"^Jvate^P*o/2(l - a;^)) - e"^~"*e-'P*o/2 (1 - a;^)"^) ' 
When Re{ip^oa^^) < 2mN we pick up the residue at e~"*^"*e*^*°/^(l— 0;^), and the amplitude 



IS 

5/ N ^ !_ -27ri(l-a;fc)~' 

2''^' Wa (1 - - e2--«*e-P- (1 _ ^,)-i 
i — 27ri 

~ 13 ^2 ^2mNate-ip*o - 1 + - Oul ' ^^^^ 

Since we are interested in imaginary p^o we switch variables, 

E + 2mN — ^p*oQ;7^, (36) 

where E is the energy in Minkowski space, with rest energy excluded. Finally we get 
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and the pole in the bubble chain sum occurs when 

^^^^ = il-6h)He- 



-Cat 



(39) 



Using ()24|1 we can determine C for various spatial and temporal lattice spacings. Results 
for the lattice spacings used in the simulations presented here are shown in Table 2. We have 
also determined which will be needed when computing the average energy by varying /? 



dat ' 



with fixed L+. 



Table 2: Contact potential coefficients (MeV 



a-i(MeV) 




C(MeV-2) 


dC 

dat 


50 


24 


-11.6 X 10-5 


-2.21 X 10-5 


60 


32 


-10.1 X 10-^ 


-2.31 X 10-5 


60 


48 


-8.75 X 10-5 


-1.91 X 10-5 


70 


64 


-7.58 X 10-5 


-1.92 X 10-5 


80 


72 


-6.96 X 10^5 


-2.04 X 10-5 



VI. BUBBLE CHAIN SUMMATION 



In this section we discuss a simple semi-analytic calculation that we use to compare with 
the results of our simulations. At T = and if kpldnnl small the energy and particle densities 
can be calculated as an expansion in kplannl- If the scattering length a„„ is small this is 
equivalent to a perturbative expansion in the coupling constant C. If a„„ is not small then 
an infinite set of particle-particle bubbles has to be summed. This is particularly obvious in 
the lattice cutoff scheme employed in this work. Since the coupling constant C is fixed by 
matching the particle-particle bubble sum to the experimental scattering length at a given 
lattice spacing, a perturbative expansion of the equation of state in powers of C will not be 
cutoff independent. An approximation scheme that will reproduce the lowest order kpann 
expansion is the bubble chain summation shown in Figs. |2 and El 

The problem at T = is that the scattering length is very large, and the expansion in 
kplann] is not useful unless the density is extremely small, p < lO-^Af- When kplannl is 
not small then corrections must be summed to all orders, and it is not obvious that there 
is any subset of diagrams that can approximate the full non-perturbative result. We note, 
however, that the bubble chain diagrams contain as a subset the diagrams with the minimum 
number of hole lines. These diagrams are summed by the low density hole line expansion. 
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oo / ' 

E-to - c 

n=0 1 2 n 



FIG. 2: Bubble chain diagrams contributing to the neutron self-energy. 



The situation is simpler if the temperature T is large compared to the degeneracy tem- 
perature Tp = (37r^p)^/^/(2m). In this case a new length scale, the thermal wavelength or 
localization length appears 



This length scale acts as an infrared regulator, cutting off long-distance correlations beyond 
this scale. In particular, it regulates the neutron-neutron scattering amplitude near thresh- 
old by giving the function B{E) in ()38|) a correction of order 0{\q^^). The net effect is that 
neutrons now have an effective scattering length of 

\aeff \ ~ min(|a„„|. At). (41) 

The expansion in a^ffP converges as long as alffP < 1 which is equivalent to T > Tp. In 
the following we compute the bubble chain diagrams shown in Figs. El and El 

The bubble chain diagrams in the neutron self-energy form a geometric series. The sum 
is given by 

where 



B 



k 

(43) 

We use this to compute the full neutron propagator 



and the average number of neutrons is 



1 - ^ ^3 ^Z}^""(fc)e-^^*° 

* k 



(45) 
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FIG. 3: Bubble chain diagrams contributing to the logarithm of the partition function. 



In a similar fashion we compute the contribution of bubble chain diagrams to the loga- 
rithm of the partition function. The relevant diagrams are shown in Fig. El The factor of 
is the due to the cyclical symmetry of the diagram. We find 

In Zg = In Zq^'^ 

1 ^ - In [l - (1 - 6hf (e"^"* - l) B{p + q, /i)] Df'-^^{p)Df'^^^{q) 

From this we can compute the average energy E 

E = -^^^ + {-niM + ^^)A, (47) 

where we have subtracted out the rest energy. The derivative with respect to [3 is calculated 
at fixed Lt by varying at, 

E = -^^^ + i-^N + f^)A. (48) 
Lt oat 

We must take into account the dependence of C on a^, and ^ for various lattice spacings 
are shown in Table 2. 



VII. COMPUTATIONAL METHODS 



We use the hybrid Monte Carlo (HMC) algorithm 2^ to generate field configurations. 
Roughly 10^ five-step HMC trajectories were run, split across 9 processors running com- 
pletely independent trajectories. Averages and errors were computed by comparing the 
results of each processor. While the HMC algorithm has become standard in lattice QCD, 
it may not be so well known in the general nuclear theory community. We therefore include 
a brief overview of the method as applied to our simulation. 
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We want to sample the partition function 

Zg(xJ DsDipDr exp [-,p;Q'..{s)i^^ - V{s)] (49) 

where the s^'s are bosonic fields and the ipi^s are fermionic fields. We use a prime since we 
will redefine Q'^j shortly. We can rewrite this as a bosonic path integral 

Zg(x J DsD(f)D(j)*ex.p[-S'{(l),s)] (50) 

where 

S'{^,s) = (P*Q^\s)<p, + Vis). (51) 

The (f)iS are bosonic fields and are called pseudofermion fields. The partition function can 
be written as 

Zgoc J DsDpD({)D(p*exp[-H{(p,s,p)], (52) 



where 



We note that 



H{(p, s,p) = S'{(f), s) + ^PaPa- (53) 



dS'{<f>,s) dQ'^'is) dV{s) 



dSa * dSa ds, 



- -<t>lQ'r,\s)^^Q'^,\s)4>i + (54) 



In our case the determinant of Q' is real and non-negative. We can therefore replace Q\j 
by a positive semi-definite Hermitian matrix Qij, with the same determinant. In our case 
Q\j has the block diagonal structure 



K 
K 



(55) 



one block for the up spins and one block for the down spins. Clearly K^j is a matrix with 
half the dimension of Q'--. If we let 

Q = K^fR (56) 

then 

Q-^ ^ K-^ {K^y' (57) 
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and 



detQ = detQ'. 
So now instead of s) we can use 



5(0,s)=0*Q-.i(s)0, + \/(s) 



where 



When we compute the derivative with respect to Sq, we have 



dsr 



Let us define 



Then 



Therefore 



,_d_ 



dK 
ds„ 



Vj - Vi 



ds„ 



dSi4>, s) 
9s„ 



dK 

dSr. 



Vj - Vi 



dSr 



+ 



dV{s) 

dSr. 



The steps for the HMC algorithm are now as follows. 



Step 1: Select an arbitrary initial real- valued configuration 

Step 2: Select a complex- valued configuration according to the Gaussian 
distribution, 

P(0)ocexp [- |e,f], 

and let 

Vi = Kr^%. 
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Step 3: Select real- valued p° according to the Gaussian random distribution 

1 



P(p°) oc exp 



Step 4: Let 



and 



for some small positive s. 



Step 5: For n = 0, 1, - 1, let 



dK 



+ 



dV{s) 

dSr, 



Sa{n + 1) = Sa{n) + Spain), 



Pa{n + l) ^Pa{n) - e 



-C 



dK 

dSr, 



ds„ 



dV{s) 
ds„ 



Step 6: Let 



Pa{N)=Pa{N) + 



dK 

dSr, 



Vj - Vi 



dK^ 

dSr, 



+ 



dVjs) 

dSry 



Step 7: Select a random number r e [0, 1). If 



< exp [-H{<l),s{N),p{N)) + H{(l),s',p')] 



then let 



s° = s{N). 

Otherwise leave s° as is. In either case go back to Step 2. 



The total number of neutrons, A, is 

1^ _ l jDsDc'Dc*f^e.p[-S] 
pdfi"^ ^ p fDsDc'Dc*exp[-S] 



2L^ 



g(m^-M)at J DsDc'Dc*c\ (H + 0)c| (n) exp [-S] ' 
/ DsDc'Dc* exp [-S] 
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for any lattice site n. Dividing by the volume V = gives the density p in lattice units. 
We compute the total energy using 

1 dlnZc 



E = + {-rriN + p)A 

Lt oat 

1 / DsDc'L'c*gexp[-5] 



Lt f DsDc'Dc*exp[-S] 
where we take into account the at dependence of C when computing We then have 



(78) 



E = 2L- 



J DsDc'Dc*f{s,c',c*)exp [-S] 
/ DsDc'Dc* exp [-S] 



+ (-mjv + p)A 



(79) 



where 



f{s, c', c*) = -(mjv - /i)e('"^-^)"*4(n + 0)c:f(n) 
d 
da. 



^\/-Cats{n)+- 



^{l-6h)]c\{n)c*{n) 



+ 



2^ X] [^(^ + ^s)c\{n) + c\{n - L)C|(n) 



(80) 



for any lattice site n. 



VIII. FREE NEUTRON RESULTS 

To better understand our lattice discretization errors, we compare our free neutron results 
on the lattice with the continuum free Fermi gas. For a continuum free Fermi gas, the 
logarithm of the partition function is 

In Zj''^^ = In Z^J + In Z^TJ' = 2 In Z^J, (81) 

where the logarithm of the single spin partition function is 



InZ!''" = V 



(27r)3 
y roo 



In 



1 + e 



/3 T5rT7+"^iV-M 



27r2 



Therefore the energy density is 

E:!'"' 1 



V 27r2 



dpp In 



dp 



(82) 



2m^ ^/3(_^+^^_/,) ^ ^ 



(83) 
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2: 

Cl 




]x-mf^ (MeV) 

FIG. 4: Density versus chemical potential for free neutrons on the lattice at T = 8 MeV and various 
lattice spacings.The curves labeled fl-f4 refer to the lattice spacings defined in Table 1. The curve 
labeled fc shows the continuum limit for free neutrons. 



and the number density is 

A free i rco -i 

We double these to get the results for both spins. In the limit as p-^**^^ we find the 
usual equipartition result for the energy per neutron, 

Tpfree o 

^ = -T. (85) 

/[free 2 ^ 

A plot of density versus chemical potential at temperature T = 8 MeV is shown in Fig. 
HI The energy per neutron at temperature T = 8 MeV is shown in Fig. El In order to avoid 
large cutoff effects, we only present results at densities corresponding with lattice fillings of 
about one-quarter or less. This is why our data at longer lattice spacings terminates at 
lower densities. 

We see in Figs. |3] and El some residual dependence on lattice spacings. While it is a small 
effect, it does make it visually confusing to overlay plots for different lattice spacings. We 
mentioned the possibility of using improved actions to reduce residual lattice discretization 
error. In this analysis, however, we use a less expensive route. We will simply rescale our 
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0.1 0.15 0.2 0.25 

P/Piv 

FIG. 5: Energy per neutron versus density for free neutrons on the lattice at T = 8 MeV and 
various lattice spacings. 



densities and energies so that the free lattice results and free continuum results agree at 
H = rriN, 

p^''''"^{a — 0,at — 0, L — oo, T, /j, = ttin) 



p(a, at, L, T, p) p{a, at, L, T, p) 



pf''''^{a,at,L,T,ijL^ uin) 



(86) 



E{a, at, L, T, p) E{a, at, L, T, p) 



E^'^^'^a = 0,at = 0,L^oo,T,p^ mjv) 



(87) 



Ef^^^{a,at,L,T,p^mN) 
We will apply the same multiphcative adjustment to all lattice results. This includes free 
lattice results, bubble chain diagram results, and lattice simulation results. 



IX. VOLUME DEPENDENCE 

For the T — 8 MeV simulations we use a lattice volume of (13 fm)^ or larger. For the 
T — A MeV simulations we use a lattice volume of (20 fm)^ or larger. The dimensions of 
our X Lt lattices are shown in Tables 3 and 4. 
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Table 3: Lattice dimensions for T = 8 MeV 





r 




1 


4 


4 


2 


4 


6 


3 


5 


8 


4 


6 


9 



Table 4: Lattice dimensions for T = 4 MeV 



Label 


L 


u 





5 


6 


1 


6 


8 



We have run simulations at both smaller and larger volumes. In Table 5 we show the results 
for T = 8 MeV and // - = -2 MeV, a'^ = 60 MeV, and a^^ = 32 MeV. Since we 
are not changing the lattice spacings for this comparison we can compare raw data without 
rescaling p and E. 



Table 5: L dependence for T = 8 MeV 



L 


pfree 


f£(MeV) 


^bubble 


5::::: (MeV) 


^simulation 


^simulation _ T 


PN 




PN 




PN 


j^sirnulation V j 


3 


0.04588 


12.582 


0.08176 


6.514 


0.0885(4) 


6.19(2) 


4 


0.04596 


12.777 


0.08215 


6.539 


0.0890(2) 


6.21(2) 


5 


0.04600 


12.757 


0.08217 


6.531 


0.0886(3) 


6.19(2) 


6 


0.04600 


12.756 


0.08217 


6.531 


0.0891(3) 


6.20(2) 



In Table 6 wc show analogous results for T = 4 MeV and ii = uin, a ^ = 50 MeV, and 
= 24 MeV. 

Table 6: L dependence for T = 4 MeV 



L 


pfree 


f£(MeV) 


pbubble 




^simulation 


^simulation ^ "\ 7"\ 


Pn 




Pn 


|™(MeV) 


Pn 


J^simulation V j 


4 


0.02234 


7.349 


0.04663 


3.474 


0.0536(3) 


3.33(2) 


5 


0.02238 


7.344 


0.04667 


3.469 


0.0533(2) 


3.33(2) 


6 


0.02238 


7.341 


0.04666 


3.469 


0.0530(2) 


3.35(2) 


7 


0.02238 


7.341 


0.04666 


3.469 


0.0530(2) 


3.35(2) 



These results suggest that finite volume effects for the lattice sizes listed in Tables 3 and 4 
are smaller than our statistical errors. If we take the volume dependence from the free and 
bubble chain calculations as a guide, then the finite volume errors are well below the 1% 
level. 
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FIG. 6: Density versus chemical potential at T = 8 MeV and various lattice spacings.The curves la- 
beled fl-f4 show free neutron results, bl-b4 show the bubble chain results, and sl-s4 show numerical 
simulations. The corresponding lattice spacings are given in Table 1. 
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FIG. 7: Density versus chemical potential at T = 4 MeV and various lattice spacings. 



X. DENSITY VERSUS CHEMICAL POTENTIAL 



In Fig. iniwe plot density versus chemical potential for T = 8 MeV, and in Fig. Owe plot 
density versus chemical potential for T = 4 MeV. In both cases we see agreement among 
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0.25 



FIG. 8: Energy per neutron versus density at T = 8 MeV for various lattice spacings. 



data for different lattice spacings. This suggests that we have properly renormalized the 
interaction and absorbed the lattice spacing dependence into the scale dependent interaction 
coefficient. We observe no phase transitions as a function of chemical potential. We note, 
in particular, that we can choose the chemical potential such that the occupation number 
in the interacting theory remains small. This implies that there is no instability towards 
a fully occupied ground state. As expected for a theory with attractive interactions the 
density at a given chemical potential is larger in the interacting theory. We observe that this 
behavior is well described by the bubble chain results for T > Tp, the low-density regime 
where we expect agreement. 

XI. ENERGY PER NEUTRON VERSUS DENSITY 

In Fig. ISlwe plot energy per neutron versus density for T = 8 MeV, and in Fig. El we 
plot energy per neutron versus density for T = 4 MeV. In both cases we see good agreement 
among data for different lattice spacings. As expected, the energy per particle approaches 
1.5T in the dilute limit. The energy per particle decreases as a function of density in the 
regime that we have studied. At small p the slope is very steep, which is consistent with 
the perturbative result for particles that have a large negative scattering length. We note, 
however, that the behavior is not linear, even at very small density. The simulations are 
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FIG. 9: Energy per neutron versus density at T = 4 MeV for various lattice spacings. 

very well described by the bubble chain results for T > Tp. We also observe that for the 
larger densities, p > O.lpAr at T = 8 MeV and p > O.OSpN at T = 4 MeV, the energy 
per particle is smaller than the result for a free neutron gas at zero temperature. In this 
regime the Fermi energy of the degenerate system is lower than the temperature. Our results 
suggest that the parameter ^ defined in (P) is smaller than 0.5. We should note, however, 
that the parameter kpTun ~ kpa-, where a is the lattice spacing, is of order 1 and simulations 
at lower temperature will be required in order to make more definitive estimates of ^. 

The decrease in the energy per neutron with increasing density does not necessarily imply 
an instability to neutron clustering. At nonzero temperature entropy must also be taken 
into account, and the question of whether or not phase separation occurs will be resolved in 
the next section when we look at the equation of state. 



XII. EQUATION OF STATE 

We integrate the density as a function of chemical potential to measure the pressure, 

p{p')dp'. (88) 
We perform the integration by least-squares fitting p{p') with a function of the form 

p{p) = (co + ci/i' + C2p'^) exp{bp). (89) 



T 1 

P = -In Zg = - A{p')dp' 
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We can then perform the integration analytically. In Fig. ^Jwe plot the pressure versus 
density for T = 8 MeV, and in Fig. ^2 we plot the pressure versus density for T = 4 MeV. 
In both cases the pressure is a smooth strictly increasing function of density. Therefore we 
conclude that there is no indication of phase separation. 

There are many models we could use to compare with our results. We first consider 
the results of a variational calculation by Friedman and Pandharipande j^. They use a 
realistic Hamiltonian that consists of the Argonne f 14 interaction supplemented by a three- 
body force. We have taken the data for different temperatures given in Table 4 of and 
interpolated to obtain the pressure for T = 4 MeV and T = 8 MeV. The result are shown 
by the crosses in Figs. El and ^2 We observe that the agreement with our calculations is 
remarkably good. There are a number of factors that are likely to contribute to this result. 
One is the fact that for the low densities considered in the present work explicit pions as 
well as three-body forces are not important. Another point is that we work with relatively 
coarse lattices. On these lattices the lattice spacing is close to the effective range parameter 
in neutron-neutron scattering. 

We also show the equation of state for a simple phenomenological model of the equation 



of state described in l23|l and the review article 



24| . The model contains a parameterization 
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P/PN 



0.12 0.14 



FIG. 11: Pressure versus density at T = 4 MeV for various lattice spacings. The crosses show the 
results of a variational calculation by Friedman and Pandharipande ^. 



of the equation of state of symmetric nuclear matter adjusted to the saturation properties 
and the compressibihty. The authors consider three possible functional forms of asymmetry, 

(version 1) Pasy = 2eaPN ( — ) (90) 

\PnJ 

(version 2) Pasy = CaPN f — ) (91) 

\PnJ 

3 

(versions) Pasy = T^^^aPN ( — ] S^, (92) 

^ \Pn/ 

where ~ 20 MeV, and the asymmetry parameter 6 is defined as 

6 = P^^^. (93) 

Pn + Pp 

The equations of state of the three different models for T = 8 MeV and T = 4 MeV are 
shown in Figs. El and El For comparison, we also show the variational results of Friedman 
and Pandharipande. We observe that our results, as well as the results of Friedman and 
Pandharipande, appear to agree most closely with version 3. Further investigations are 
needed to determine whether other properties of this simple model agree with our lattice 
simulation results. 
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FIG. 12: Pressure versus density at T = 8 MeV for the phenomenological equation of state discussed 
in [2^. The three different curves correspond to different parametrizations of the symmetry energy. 



XIII. SUMMARY AND CONCLUSIONS 



In this work we studied neutron matter by combining pionless effective field theory at 
lowest order with non-perturbative lattice methods. To determine the neutron contact 
interaction we summed bubble chain diagrams contributing to neutron-neutron scattering 
at a given lattice spacing. The contact interaction was then adjusted to produce the pole in 
the amplitude indicated by Liischer's finite volume formula for the physical ^Sq scattering 
length. Having determined the interaction coefficient for various lattice spacings, we then 
simulated neutron matter on the lattice using hybrid Monte Carlo at temperatures 4 and 8 
MeV and densities below one-fifth normal nuclear matter density. 

We find that our results at different lattice spacings agree with one another. This suggests 
that the continuum limit exists and that our effective theory was properly renormalized or, 
more conservatively, that any cutoff dependence is numerically small. For the range of 
parameters studied in this work we observe no instabilities towards phase separation, or 
towards lattice artifacts such as a completely filled lattice. While not unexpected, this is 
not a trivial result since the non-perturbative simulation includes all possible diagrams. 

The energy per particle at temperatures T = 4 MeV and T = 8 MeV shows a steep 
downward slope at very small density, which is a sign of the strong attractive interaction 
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FIG. 13: Pressure versus density at T = 4 MeV for the phenomenological equation of state discussed 
in [2^. The three different curves correspond to different parametrizations of the symmetry energy. 



between neutrons. At intermediate densities p ~ 0-lpN the energy per particle levels off. This 
behavior is reproduced quantitatively by our bubble chain calculations for T > Tp. In the 
future we wish to push our simulations to lower temperatures and determine the universal 
parameters ^ and C for the energy per particle and gap as defined in (P). Simulations in this 
regime will require a source term for the di-neutron field. We have not seen unambiguous 
signs of superfluidity in our simulations. We did observe, however, a significant drop in 
HMC acceptance rate for the simulation at the lowest temperature and highest density 
studied in this work. This may well be an indication for the onset of superfluidity. 

Our results for the pressure of pure neutron matter agree remarkably well with the vari- 
ational calculation of Friedman and Pandharipande j^. In the future it will be interesting 
to study whether the agreement persists if higher order terms in the effective Lagrangian or 
explicit pions are introduced. We have also performed simulations with explicit pions jl7|. 
but these data were taken at larger density and temperature. It will also be interesting to 
study systems with a finite proton fraction. This is easiest in the limit of exact Wigner 
symmetry, as the leading order Euclidean action is positive in that case 25 1. 

It will also be interesting to investigate the phase structure of the leading order effective 
theory in more detail. In order to have a positive Euclidean action the coefficient C of the 
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four-fermion interaction has to be negative. This corresponds to either a negative scattering 
length or scattering length that is large and positive [l8|. This implies that the effective 
theory studied in this work can be used to investigate the BCS-BEC crossover in a dilute 
Fermi gas. 
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